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We apply the atom-atom potentials to molecular crystals of iron (II) complexes with bulky 
organic ligands. The crystals under study are formed by low-spin or high-spin molecules 
of Fe(phen)2(NCS)2 (phen = 1,10-phenanthroline), Fe(btz)2(NCS)2 (btz = 5,5',6,6'-tetrahydro- 
4_f/,4'_ff-2,2'-bi-l,3-thiazine), and Fe(bpz)2(bipy) (bpz = dihydrobis(l-pyrazolil)borate, and bipy = 
2,2'-bipyridine). All molecular geometries are taken from the X-ray experimental data and assumed 
to be frozen. The unit cell dimensions and angles, positions of the centers of masses of molecules, 
and the orientations of molecules corresponding to the minimum energy at 1 atm and 1 GPa are 
calculated. The optimized crystal structures are in a good agreement with the experimental data. 
Sources of the residual discrepancies between the calculated and experimental structures are dis- 
cussed. The intermolecular contributions to the enthalpy of the spin transitions are found to be 
comparable with its total experimental values. It demonstrates that the method of atom-atom 
potentials is very useful for modeling organometalic crystals undergoing the spin transitions. 



I. INTRODUCTION 



The Crystal Field Theory (CFT), proposed in [Ij] and 
known to majority of chemists through suggests that 
coordination compounds of d-elements with electronic 
configurations d^, (f" , cP or (f can exist either in high- 
spin (HS) or low spin (LS) forms (sometimes interme- 
diate values of the total spin are also possible). In the 
case of strong-field ligands the d-level splitting measured 
by the average crystal field parameter IQDq exceeds the 
average Coulomb interaction energy of d-electrons P and 
the ground state is LS. In the case of weak-field ligands 
with IQDq <^ P, the ground state is bound to be HS. 
If, however, lODq = P, the LS and HS forms of the 
complex may coexist in equilibrium, and the fraction 
of either spin form depends on temperature, pressure, 
and/or other macroscopic thermodynamic parameters. 
The process when the fraction of molecules of different 
total spin changes due to external conditions is called a 
spin crossover (SC) transition. For the first time this 
phenomenon was reported in 1931 3. Nevertheless, ex- 
tensive studies of SC started only in 1960s- 70s. Nowa- 
days, dozens of complexes capable to undergo spin tran- 
sitions (spin-active complexes) are known, and most of 
them are those of Fe(II). A general review of the field 
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can be found in [J| . 

A wealth of potential practical applications like dis- 
lays and data storage devices (see a detailed review in 
) is one of the reasons for research activity in this area. 
Industrial applications pose strict demands on the char- 
acteristics of the materials to be used. As a consequence, 
the problem of predicting SC transition characteristics 
(whether it is smooth or abrupt, the transition temper- 
ature, the width of the hysteresis loop, the influence of 
additives [^) is of paramount importance. Theoretical 
description of spin transitions is a great challenge by it- 
self, and until now a coherent theory allowing to relate 
the composition of the materials with the characteristics 
of the transition has not been developed. Discussion of 
these issues and an overview of the existing theories are 
given in Q. 

In general, the SC modehng includes two aspects: (i) 
that of the interactions within one molecule of a spin- 
active complex, and (ii) that of the interactions between 
these molecules. The latter is crucially important for un- 
derstanding of specific features of the SC transitions in 
solids because the SC manifesting itself as a first-order 
phase transition is controlled by intermolecular interac- 
tions. These ideas are built in the simplest model ca- 
pable of describing spin transitions in solids proposed 
by Slichter and Drickamer [8|. This model considers the 
solid as a regular solution of molecules in the LS and HS 
states. The model predicts, in agreement with the exper- 
iments, that the spin transition may be either smooth or 
abrupt or may exhibit hysteresis, and its character is de- 
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termined by a phenomenological intermolecular parame- 
ter r, specific for each materiaL However, the experimen- 
tal data on the heat capacity and the X-ray diffraction 
contradict to this model. 

The thermal dependence of the heat capacity of the 
Fe(phen)2(NCS)2 crystal is better explained by an alter- 
native domain model Q. Diffraction patterns of spin 
transition crystals, measured at intermediate temper- 
atures, simultaneously contain the Bragg peaks corre- 
sponding to the pure LS and HS phases, while no peaks 
for intermediate lattice of a solution were observed p^ . 
Another problem is that the parameter F is phenomeno- 
logical one, and it cannot be sequentially derived in 
terms of microscopic characteristics of the constituent 
molecules or their interactions. At the same time, within 
the Slichter-Drickamer model, the type of behavior is 
tightly related to the sign and magnitude of F, so that 
a smooth transition requires F > 0, an abrupt transi- 
tion occurs at F < and hysteresis is possible only if 
F < is less than some critical threshold, which in its 
turn depends on the transition temperature [Sj. It has 
been shown that if the relaxation of the lattice is not al- 
lowed, then under very natural assumptions F is positive 
11 1 , but the lattice relaxation can lead to F of either sign 

Significant progress in the understanding of the spin 
transitions in crystals is attributed to the Ising-like mod- 
els of intermolecular interactions in spin-active materials 
pljj . Adaptations of the initial Ising model to the spin 
transitions include corrections for intramolecular vibra- 
tions, domain formation, parameters distribution, elastic 
distortions, presence of two metal atoms in a spin-active 
molecule, etc. 0, [13] ■ These models do not have an- 
alytical solutions and they are solved either in a mean 
field approximation which leads to results analogous to 
(or even coinciding with) the Slichter-Drickamer model 
\j\ or numerically. 

In spite of the diversity of the models used in the liter- 
ature, the theoretical description of the spin transitions 
is not yet satisfactory. First, the existing theories are not 
capable to reproduce the whole set of the experimental 
data (e.g. asymmetry of the hysteresis loop Q). Second, 
all of them contain phenomenological parameters, like F 
in the Slichter-Drickamer model, or the energy gap Aj 
or the interspin interaction constants Jij in the Ising-like 
models, or the bulk modulus K and the Poisson ratio a 
in (l5| {K and a can be measured, but for the purpose 
of the theory they must be independently predicted), 
etc. Third, even if the models include microscopic level 
consideration, they use oversimplified description of the 
molecules (as spheres, ellipsoids), which is not sufficient 
for constructing a complete theory, especially due to im- 
portance of the short intermolecular contacts tentatively 
responsible for the cooperativity effects (tt-tt interactions, 
S- ■ • H— C interactions, hydrogen bonds, etc. jl6|). 

These shortcomings can be overcome by using explicit 
potentials for intra- and intermolecular interactions. In 
this case one may expect to obtain independent estimates 



of the numerous parameters required by the phenomeno- 
logical theories. These potentials should also be a help- 
ful tool for checking the validity of the initial postulates, 
such as the formation of a regular solution or the domain 
structure, thus clarifying some obscure points in the the- 
ory itself. 

An adequate ah initio calculation of the energies of iso- 
lated transition metal complexes, and moreover those of 
the crystals formed by these complexes, is a very compli- 
cated problem. Significant electron correlation within the 
d-shells breaks the self-consistent field approximation, so 
that explicit account of nontrivial (static) electron corre- 
lation is unavoidable. The existing implementations of ah 
initio approaches for solids fail to provide the necessary 
quality of the results. 

There is a number of attempts to use the DFT-based 
methods to take into account the electron correlation in 
the SC complexes [l3| . These methods yield good results 
for many characteristics of isolated spin-active molecules 
(optimal molecular geometry, Mossbauer parameters, vi- 
brational frequencies, nuclear inelastic scattering spec- 
tra) [ll,|Ti,|2fl,|2l|. However, the DFT in its traditional 
form, as it is demonstrated in is not capable to re- 
produce coherently the static correlations, which are ex- 
tremely important for the correct description of the spin 
transitions even in an isolated molecule. For that rea- 
son the results for the energy gap between the LS and 
HS states, and hence for the transition temperature, ob- 
tained by the DFT techniques are absolutely disastrous. 
The common versions of DFT, such as B3LYP, often pre- 
dict a wrong ground state multiplicity, let alone the value 
of the energy difference fl7'|. For example, the tempera- 
ture of the spin transition in Fe(phen)2(NCS)2 was found 
to be an order of magnitude too large (1530 K instead 
of 176 K) ^23;]. In addition, most DFT studies are lim- 
ited to isolated molecules in vacuo, and the heat of the 
spin transitions in a crystal is identified with the energy 
difference of isolated molecules. The influence of the in- 
termolecular interactions is thus neglected. 

Only a few isolated attempts to explicitly model 
a spin-active crystal by the DFT method have been 
reported 0, 111, IH [Hi. The application of the 
LDA approximation with the periodic boundary con- 
ditions to the crystals of Fe(trim)2X2 (X = F, CI, 
Br, or I, and trim = 4-(4-Imidazolylmethyl)-2-(2- 
imidazolylmcthyl)imidazole) formed by either LS or HS 
molecules (20,] demonstrated that the intermolecular in- 
teractions strongly affect the energy splitting between 
the LS and HS isomers, thus necessitating their adequate 
treatment within coherent SC models. The experimen- 
tal X-ray structures for some complexes are available, so 
that the optimal geometry of the crystals found by LDA 
can be verified. This comparison showed that the unit 
cell volumes were overestimated by 20-24%. At the same 
time, the calculated N...X distances were 0.1-0.3 A lower 
than the experimental ones, and the 7r-stacking distances 
were underestimated by 0.3-0.7 A. Although in general 
it is difficult to separate the errors from intra- and inter- 



3 



molecular interactions, the geometry of individual spin 
isomers is usually described much better than the rela- 
tive position of the molecules in the crystal. The GGA 
approximation has been used to optimize molecular ge- 
ometry and lattice parameters of the LS and HS crys- 
tals of [Fe(pyim)2(bipy)l(C104)2-2C2H50H (pyim = 2- 
(2-pyridyl)imidazole) |26j . The bond lengths were found 
to be quite reasonable. However, the lattice parameters 
were poorly reproduced, so that even the wrong sign of 
the unit cell volume change for the LS to HS transition 
was obtained: (—6.82 A"^ instead of experimental value 
of -h228.02 A^). The authors explain it by "the well- 
known shortcomings of DFT methods in application to 
weak intermolecular interactions" [1^. The DFT+U ap- 
proach with the GGA approximation has been applied 
to model spin-active crystals of Fe(phen)2(NCS)2 and 
Fe(btr)2(NCS)2(H20) (btr = 4,4'-bis-l,2,4-triazole) [H. 
The studies of the Fe(phen)2(NCS)2 crystal have demon- 
strated that DFT is capable of reproducing the lattice 
parameters with the precision of 1-5% and the unit cell 
volume with the precision of 2 — 7% [131 . Unfortunately, 
these works do not employ a much better substantiated 
approach to the DFT-based treatment of van der Waals 
interactions previously proposed by the same authors, 
based on explicit treatment of correlations coming from 
the long range part of the electron-electron interactions 

Summarizing, DFT models either produce poor re- 
sults for spin-active complexes or require parameters (like 
DFT+U) adjusted to reproduce the experimental data. 
At the same time, the very idea of modeling such com- 
plex system as a crystal formed by spin-active transition 
metal complexes at a uniform level of theory seems to 
be incorrect. The systems under consideration consist 
of numerous components, and it is much more natural 
to treat these components separately - each at the ade- 
quate level of the theory. The most important separation 
is that on intra- and intermolecular interactions. On the 
level of molecules one can further separate a highly cor- 
related d-shell from the rest of the molecule. This idea 
has been implemented as a specialized quantum chemical 
method - Effective Hamiltonian of Crystal Field (EHCF) 
(29*1 which has been successfully applied to describe the 
spin isomers of Fe(phen)2(NCS)2 ^SOJ- Furthermore, it 
has been demonstrated that the geometry of spin-active 
complexes can be adequately described by the EHCF 
technique with ligands treated by molecular mechanics 
force fields [Slj . 

On the level of interactions between molecules the 
paramount fact is that the molecular crystals formed by 
spin-active molecules consist of complexes with bulky or- 
ganic ligands. Intermolecular contacts in such crystals 
are those between the organogenic atoms like C, H, N, 
S, etc. The d-shells of the central ions are effectively 
shielded by the ligands. Thus, it is reasonable to assume 
[m that the d-shells do not directly affect the interac- 
tions between the molecules of the different total spin in 
the crystal, but influence it indirectly: through the vari- 



ation of the equilibrium interatomic distances Fe— N in 
these complexes, which is further translated into differ- 
ent "sizes" of the LS and HS isomers. In this context, 
the standard methods developed for organic molecular 
crystals can be successfully applied in this case as well. 
The main purpose of the present work is to identify an 
adequate way to model intermolecular interactions for 
crystals formed by spin-active molecules. 



II. ATOM-ATOM POTENTIALS METHOD 

In order to avoid unnecessary complications, we limit 
our task in the present paper to checking the possibil- 
ity of applying the simplest method of modeling inter- 
molecular interactions - atom-atom potentials [32l | - to 
crystals formed by spin-active complexes. The method 
assumes that the energy of the molecular crystal (calcu- 
lated relative to the system of isolated molecules) can be 
represented as: 
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where each term is the energy of the interaction between 
the a-th atom of the m-th molecule in the unit cell num- 
ber r = {ra,riy,rc) and the a'-th atom of the m'-th 
molecule in the unit cell r' depending on the distance 
R. Due to the equivalence of all unit cells, we can get rid 
of summation over r', and the energy per molecule u can 
be written as: 
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where M is the number of molecules per unit cell. 

A number of approximations have been suggested for 
the atom-atom interaction. The most widespread ones 
are the Buckingham potential {&-exp): 
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and the Lennard- Jones potential (6-n) 

Eaa'(R) = 
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In the above formulae the Aaa' and Baa' parameters 
for the interaction between atoms of different types are 
often calculated as the geometric mean values of the cor- 
responding homogeneous interaction parameters: 



-^aa' \l -^aa-^a' a' , Baa' — \j BaaBa' a' , (5) 

while the Caa' parameter is approximated in a similar 
way as an arithmetic mean value: 
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Due to these approximations, the energy can be rep- 
resented as a fast computable function depending on the 
lattice parameters and relative positions and orientations 
of the molecules in the unit cell provided that molecular 
geometry of the complex is fixed. Having found the min- 
imum of this function, one gets estimates of the inter- 
molecular interaction energy (sublimation energy), the 
equilibrium unit cell parameters, and the positions and 
orientation of the molecules in the unit cell at the abso- 
lute zero temperature and absence of external pressure. 

One can easily extend the method to account for the 
external pressure. For this purpose one should optimize 
the enthalpy H instead of the potential energy U. The 
enthalpy is defined as 

H = U + PV, (7) 

where the volume V is determined by the lattice parame- 
ters. As for the thermal dependence of the lattice param- 
eters, the matter is not so simple. One should basically 
minimize the Gibbs energy G to estimate the equilibrium 
values of the lattice parameters at a non-zero tempera- 
ture (and pressure). This procedure includes calculation 
of the entropy of the crystal undergoing the spin transi- 
tioii, which is a separate no n- trivial challenge, as shown 
in [35]. To avoid this, one may confine to minimiza- 
tion of the internal energy U or the enthalpy H , but the 
resulting lattice parameters will be relevant only for the 
absolute zero of temperature. On the other hand, in prac- 
tice the parameters of atom-atom interaction are fitted 
in such a way that the lattice parameters corresponding 
to the minimum of the model internal energy U best re- 
produce the experimental lattice structures measured at 
the room temperature (see e.g. [36l|). In this case the 
model includes the entropy factor implicitly, and the lat- 
tice parameters found by direct minimization of U should 
actually refer to the room temperature. 

The accuracy of the atom-atom approach is corrobo- 
rated by extensive statistics obtained for organic molec- 
ular crystals [H, [H, [s^l . Typically it provides the ac- 
curacy level of ca. 0.1-^5 kcal/mol in energy terms for 
a wide range of organic crystals. However, in theory we 
can expect much better precision for the relative energies 
of the crystals undergoing the spin transition, since the 
LS and HS crystals are very similar to each other (as is 
shown below, the shortest contacts are the same). 



III. MODELING METHOD 

We performed calculations for the molecular crystals 
formed by each of the spin isomers of Fe(phen)2(NCS)2, 
Fe(btz)2(NCS)2, and Fe(bpz)2(bipy). The hgands are de- 
picted in Fig. [1] and the molecules themselves are shown 
in Figs. [2]|4l The objects were chosen based on the fol- 
lowing considerations. First, all these crystals consist of 
neutral molecules only, without ions or solvents. As a re- 
sult, the molecules are held together in the crystal by the 
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FIG. 1: Structure formulae for the ligands of the spin-active 
complexes studied. 




FIG. 2: Molecular structure of Fe(phen)2(NCS)2. 



van der Waals forces (no strong Coulomb forces or ob- 
vious hydrogen bonds are involved), which dramatically 
simplifies modeling of the energy. Second, these three 
substances represent all main types of spin transitions: 
abrupt one in the Fe(phen)2(NCS)2 crystal, smooth one 
in the Fe(btz)2(NCS)2 crystal, and the transition with 
hysteresis in the Fe(bpz)2(bipy) crystal. Finally, the crys- 
tallographic data (including the molecular geometries) 
for both HS and LS forms of these three substances are 
available in the literature. 

The energy of van der Waals interactions was described 
by the Lennard-Jones (6-12) potential with the parame- 
ters of the " Universal Force Field" (UFF) parameteriza- 
tion [s?! and by the Buckingham {Q-exp) potential with 
the parameters provided in [36| (see Tables HI |TT| . In the 
latter case the parameters for the C- ■ • H, N- ■ • H, S- ■ • H, 
C- • • N, and S- • ■ C interactions are given in [36j explicitly 
and there is no need to use eqs. (O and Unfortu- 
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FIG. 3: Molecular structure of Fe(btz)2(NCS)2. 




FIG. 4: Molecular structure of Fe(bpz)2(bipy). 

nately, the system of parameters [sB] for the {6-exp) po- 
tential has not been extended to boron. Hence we took 
the minimum depth and the interatomic separation at 
the minimum for the B- ■ • B pair from [3^1 , estimated the 
corresponding A, B and C parameters and found the pa- 
rameters for the B- • • H, B- • • C, B- • • N, and B- • • S pairs 
following eqs. ([5]) and ([6]). The parameters of the {6-exp) 
potential for pairs involving Fe atom(s) are not deter- 
mined, but they arc immaterial in the present context, 
and we set them to be equal to zero. 

The MOLCRYST program suite ^ capable of calcu- 



TABLE I: Parameters of the Lennard- Jones (6-12) potential 
[S^ used in the calculations. 





H 


B 


C 


N 


S 


Fe 


A, kcal-A'^/mol 


50.9 


1668 


685 


332 


2365 


15.9 


B, lO'-kcal-A'Vmol 


0.147 


38.6 


11.2 


3.99 


51 


0.048 



lation and minimization of molecular crystals energy and 
enthalpy with use of the Lennard-Jones and Bucking- 
ham atom-atom potentials was employed. This program 
has been thoroughly tested on the examples of molecu- 
lar crystals of aromatic hydrocarbons. The geometries 
of HS and LS forms of the complexes were taken from 
experiments [iol . l4ll . and were assumed to be fixed 
(frozen) throughout the modeling. The validity of the 
rigid-body approximation can be tested [l^j and the anal- 
ysis of the difference vibrational parameters for an SC 
crystal demonstrated [S] that the non-rigidity is rela- 
tively small for both HS and LS forms. 

When calculating the energy according to eq. we 
restricted ourselves to summation over three layers of 
unit cells around the central "0-th" unit cell. In other 
words, only those r = (?'a,''6,fc) were included into the 
sum, for which \ra\ < 3, \rb\ < 3, and jrd < 3. It was 
found that extending this limit to 4 or more layers does 
not affect the final result for the energy or enthalpy (the 
differences are less than 0.01 kcal/mol). As for the equi- 
librium values of the lattice parameters, their values are 
stable (within 0.1%) already with one layer of the sur- 
rounding unit cells (those adjacent to the "0-th" cell) 
included into the summation. 

To find the equilibrium values of the lattice parame- 
ters, positions of the centers of masses (CM), and the 
rotation angles of molecules in the unit cell, minimiza- 
tion of the enthalpies of six pure crystals (three HS 
and three LS) was performed. Pressure was set to be 
1 atm. In all the cases the experimental crystallographic 
data were taken as initial approximations. At the first 
stage we minimized the enthalpy as a function of five or 
six parameters (a, 6, c, one non-trivial rotation angle, 
and one non-trivial CM coordinate in the cases of the 
Fe(phen)2(NCS)2 and Fe(btz)2(NCS)2 crystals; the same 
plus the unit cell angle /3 in the case of the Fe(bpz)2(bipy) 
crystals) , preserving the symmetry of the crystal (Pbcn, 
Pbcn and C2/c correspondingly); after that we checked 
that the final point of the previous step is the global min- 
imum, allowing for variation of all 27 parameters (a, b, 
c, three unit cell angles, three rotation angles for each of 
four molecules in the unit cell, three CM position coordi- 
nates for three out of four molecules in the unit cell; the 
fourth molecule position is not independent due to the 
crystal translational symmetry). The optimized struc- 
tures are shown on Figs. [5][7l 

The enthalpy was calculated as the sum of the internal 
energy and the product of the pressure and the volume of 
the crystal with 1 mole of molecules. In all the cases the 
internal energies were found to be about —50 kcal/mol 
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TABLE 11: Parameters of the Buckingham {6-exp) potential [s^ used in the calculations. 





H-- - H 


C •H 


N •H 


S---H 


C---C 


C- ■ - N 


S---C 


N---N 


S---S 


B- • ■ B [38] 


A, kcal-A"/mol 


26.1 


113 


120 


279 


578 


667 


1504 


691 


2571 


3.688 


B, 10^-kcal/mol 


5.774 


28.87 


54.56 


64.19 


54.05 


117.47 


126.46 


87.3 


259.96 


19.84 


c, A-^ 


4.01 


4.10 


4.52 


4.03 


3.47 


3.86 


3.41 


3.65 


3.52 


6.82 



Note: parameters for the S- • ■ N interaction and interactions of 
B with other atoms were calculated according to the superposition 
approximation JSj and 




FIG. 5: Crystal structure of Fe(phen)2(NCS)2. 




FIG. 6: Crystal structure of Fe(btz)2(NCS)2. 




FIG. 7: Crystal structure of Fe(bpz)2(bipy). 



relative to the isolated molecules. The experimental data 
to verify this result are not available. However, the en- 
ergy magnitude is quite reasonable in comparison with 
the available data on organic molecular crystals [H, , 
taking into consideration that the numbers of interatomic 
contacts per molecule in the crystals under study are a 
few times higher than those in ordinary organic crystals. 
The differences between the internal energies and the en- 
thalpies in all the cases at 1 atm are rather small, less 
than 0.01 kcal/mol, which is not surprising, since we deal 
with solid substances. 

As mentioned above, the {6-exp) potential parameter- 
ization from [36j implicitly includes the entropy contri- 
bution since it was fitted to reproduce the room temper- 
ature geometries of crystals by minimization of the inter- 
nal energy rather than the Gibbs energy. So the results of 
our calculation with the Buckingham potential should be 
compared with the room temperature experimental data. 
The matter is not so clear in the case of the UFF parame- 
ter system [s^l- The authors introduce their parameters 
of the van der Waals interaction explicitly referring to 
ionization potentials, polarizabilities, and Hartree-Fock 
calculations, so these parameters seem to be providing 
unadjusted estimates of the internal energy. Neverthe- 
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less, direct comparisons of the numbers produced with 
their empirical parameters and experimental geometries 
are widely used. Strictly speaking, we do not have suffi- 
cient information to judge whether our results obtained 
with this (6-12) potential describe physical properties for 
the absolute zero temperature or for the room tempera- 
ture. However, comparing the experimental data on the 
lattice parameters of Fe(phen)2(NCS)2 at 15 K, 32 K, 130 
K and 298 K [i^ with the results of our calculations, we 
can see that the latter are somewhat closer to the high- 
temperature values of the lattice parameters, rather than 
to the low-temperaturc ones. 

The room temperature crystallographic data are avail- 
able only for the HS crystals. As for the LS crystals, 
we need to extrapolate their experimental lattice param- 
eters to the room temperatures to make the comparison 
with the results of our calculations possible. This is es- 
pecially important for the analysis of the changes of the 
lattice parameters Ay, Aa, A6, Ac, etc. in the course 
of the spin transition, otherwise the calculated experi- 
mental values would include not only the contribution 
of the spin transition itself, but also of thermal expan- 
sion of the crystal. In the cases of the Fe(phen)2(NCS)2 
and Fe(btz)2(NCS)2 compounds, the dependences of the 
V, a, b, c parameters and the HS molecules fraction x 
(from the magnetic susceptibility data) on temperature 
are known in the range from ca. 130 K to 293 K (4l| 
(each series consists of 22-25 observations). In a linear 
approximation. 



V (T, x{T)) ^{l~x) iVo,LS + ^ivMT - To)) 

+X {Vo^HS + KV,HS{T -To)) , 



and similarly for the a, b and c parameters. We deter- 
mined the coefficients Voxs-, kv,ls, VoMS, i^v,hs, ^tc. by 
the method of least squares [E? of such models are typi- 
cally 0.995-^0.999), and made extrapolation of the lattice 
parameters to the room temperature and unchanged frac- 
tion of the HS molecules. These extrapolated values were 
used for comparison with the results of the method of 
atom-atom potentials. As for the Fe(bpz)2(bipy) crystal, 
the lattice parameters, published in the literature, were 
measured only at few temperatures [13, Thermal 
coefficients of expansion, calculated for the LS form on 
different temperature intervals, differ significantly, which 
does not allow for a reliable extrapolation of the lattice 
parameters to the room temperature. On the other hand, 
high-temperature coefficients of expansion are more sta- 
ble. Because of these reasons, we extrapolated the HS 
crystal lattice parameters to 139 K to estimate AF, Aa, 
A6 and Ac free of thermal distortions, though only at 
139 K. The results of the extrapolations made are used 
in the next Section for comparison with calculated opti- 
mal lattice parameters. 



IV. RESULTS AND DISCUSSION 
A. Crystal geometries 

The most important for thermodynamical description 
of the spin transitions characteristic of the lattice is the 
unit cell volume. The estimates of this quantity obtained 
by the atom-atom potentials model are given in Tables 
IIIHTVl The average error in the computed volume is 1.8%, 
ranging from 0.5% to 4.0%. The Lennard-Jones and the 
Buckingham potentials provide comparable levels of ac- 
curacy. These numbers should be compared with the 
discrepancy of 20 - 24% in ^23] and 1 - 8% in ^26^ (both 
calculated with the DFT method), the only analogues 
published so far. At the same time one should remember 
that these data include relatively small errors in the ge- 
ometries of separate molecules while our calculations are 
free of them because we used the experimental structures 
for the molecules. 

The changes of the unit cell volumes in the course of 
the spin transition are relatively small differences of two 
large numbers, and their correct estimation is difficult. 
For example, AV of [Fe(mam)2(bpy)] (€104)2 •2C2H5OH 
was found to be negative [26] , though all complexes stud- 
ied experimentally have positive Ay, in agreement with 
the fact that the Fc— N bonds are longer in the HS com- 
plexes, and thus the HS molecules should have a larger 
"size". The calculated value of Ay for [Fe(trim)2]Cl2, 
published in [20| . has the correct sign, but the experi- 
mental volumes of the LS and HS crystals are available 
only for different temperatures, which makes it impossi- 
ble to compare the experimental and calculated values. 

The values of Ay of the Fe(phen)2(NCS)2 compound, 
calculated by us with both Lennard-Jones and Bucking- 
ham atom-atom potentials, are fairly close to the ex- 
perimental values (extrapolated to the room tempera- 
ture), being probably overestimated by 10 — 17% (while 
the uncertainty in the extrapolated experimental value 
is ca. 10%). In the case of Fe(btz)2(NCS)2, the er- 
rors are correspondingly about -1-15% and —3% for the 
two potentials, while the uncertainty in the extrapolated 
experimental value is ca. 7%. Finally, in the case of 
Fe(bpz)2(bipy) the calculated values differ from the ex- 
perimental one, extrapolated to 139 K, by 4 — 14%. We 
would like to stress that the temperature dependence of 
Ay is much stronger, than that of the unit cell vol- 
ume V. For example, the low-temperature (at 15 K) 
Ay(Fe(phen)2(NCS)2) equals to 61.2 A^, the Ay value 
extrapolated to 293 K is about 70-^79 A'^, and the differ- 
ence between the experimental unit cell volume of the HS 
form at 293 K and that of the LS form at 130 K is 119.1 
A'^. The presumable errors of the atom- atom potentials 
method in calculations of Ay (ca. 5 — 15 A'^) are compa- 
rable with the uncertainties in the extrapolated estimates 
for experimental values (ca. 3 — 9 A'^) and much less than 
the changes in the volumes of the crystals caused by tem- 
perature expansion of the crystals (dozens of A'^). 

As for the unit cells themselves, in all the cases the 
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TABLE III: Comparison of experimental and calculated unit cell parameters for Fe(phen)2(NCS)2 (at 1 atm). 





a, A 


b, A 


c, A 




V, 


angle," 


CM y/b 


H, kcal/mol 


The LS isomer 


calc. (6-12) 


13.185 


9.922 


17.347 


90 


2269.2 


142.85 


0.1112 


-54.46 


calc. {6-exp) 


12.992 


9.861 


17.281 


90 


2214.0 


144.48 


0.1138 


-54.38 


calc. {6-exp) modif. 


13.017 


9.991 


17.469 


90 


2271.7 


144.59 


0.1065 


-54.04 


exp. 15 K [46] 


12.762 


10.024 


17.090 


90 


2186.3 


143.84 


0.0943 




exp. 130 K [40] 


12.770 


10.090 


17.222 


90 


2219.1 


140.51 


0.0925 




exp. extrap. to 293 K 


12.77 


10.18 


17.40 


90 


2259 








The HS isomer 


calc. (6-12) 


13.525 


9.910 


17.583 


90 


2356.7 


147.36 


0.1071 


-52.65 


calc. {6-exp) 


13.264 


9.869 


17.542 


90 


2296.2 


149.35 


0.1088 


-53.80 


calc. {6-exp) modif. 


13.227 


10.017 


17.815 


90 


2360.4 


149.19 


0.0993 


-52.50 


exp. 15 K [46] 


13.185 


9.948 


17.135 


90 


2247.5 


153.84 


0.0989 




exp. 293 K [40] 


13.161 


10.163 


17.481 


90 


2338.2 


147.09 


0.0938 




The difference between the HS and LS isomers 


calc. (6-12) 


0.340 


-0.012 


0.236 





87.5 


4.51 


-0.0040 


1.81 


calc. {6-exp) 


0.272 


0.008 


0.260 





82.2 


4.87 


-0.0050 


0.57 


calc. {6-exp) modif. 


0.210 


0.026 


0.347 





88.7 


4.60 


-0.0072 


1.54 


exp. extrap. to 293 K 


0.39 


-(0.024-0.04) 


0.05H-0.08 





70-^79 








exp. (15 K) [46| 


0.423 


-0.076 


0.045 





61.2 


10.00 


0.0045 





Note: see detailed explanation of "{6-exp) modif." parameteri- 
zation in Subsection IIIV Bl l. 



TABLE IV: Comparison of experimental and calculated unit cell parameters for Fe(btz)2(NCS)2 (at 1 atm). 





a, A 


b, A 


c, A 


p; 


V, 


angle," 


CM y/b 


H, kcal/mol 


The LS isomer 


calc. (6-12) 


13.266 


10.518 


16.889 


90 


2356.4 


125.60 


0.0385 


-54.15 


calc. {6-exp) 


13.099 


10.498 


16.741 


90 


2302.1 


126.45 


0.0445 


-57.10 


calc. {6-exp) modif. 


13.160 


10.652 


16.963 


90 


2377.9 


127.77 


0.0493 


-52.55 


exp. (130 K) [41] 


13.055 


10.650 


16.672 


90 


2318.1 


127.48 


0.0421 




exp. extrap. to 293 K 


13.17 


10.80 


16.88 


90 


2397 








The HS isomer 


calc. (6-12) 


13.242 


10.724 


16.947 


90 


2406.6 


129.45 


0.0451 


-54.32 


calc. (6-exp) 


13.077 


10.669 


16.803 


90 


2344.3 


130.14 


0.0498 


-58.20 


calc. (6-exp) modif 


13.190 


10.786 


16.973 


90 


2414.5 


130.77 


0.0527 


-53.61 


exp. (293 K) [41] 


13.288 


10.861 


16.920 


90 


2441.9 


129.79 


0.04150 




The difference between the HS and LS isomers 


calc. (6-12) 


-0.023 


0.206 


0.059 





50.2 


3.85 


0.0066 


-0.17 


calc. (6-exp) 


-0.022 


0.172 


0.062 





42.2 


3.68 


0.00535 


-1.10 


calc. (6-exp) modif. 


0.030 


0.134 


0.010 





36.6 


3.00 


0.0034 


-1.06 


exp. extrap. to 293 K 


0.12 


0.06 


0.03-^0.04 





42^45 









Note: sec detailed explanation of "{6-exp) modif." parameteri- 
zation in Subsection HIV Bl l. 



symmetry for the energy minimum points, according to 
our calculations, coincides with the experimental one. 
Orientation of a molecule in the unit cell can be charac- 
terized by three angles, corresponding to the transforma- 
tion of coordinates from the molecular coordinate system 
{e.g. that of the principal axes of inertia tensor) to the 
laboratory (or crystal) coordinate system. In all the con- 
sidered cases, two of these angles have trivial values (0, 
90 or 180°); the values of the third angle, corresponding 
to rotation around the C2 axis of the molecule, are given 
in Tables IIIHTVl The same is true for the CM positions 



of molecules within a unit cell. Two parameters out of 
three for each molecule are trivial (0, 1/4, 1/2, or 3/4 
of the corresponding translation period). The remaining 
parameter (corresponding to the y coordinate in the units 
of b) is given in Tables IllUfVl as well. One can see that 
the calculated values are fairly close to the experimental 
ones both for the rotation angles and the CM positions. 

The discrepancy between the calculated and experi- 
mental values of the lattice parameters a, b, c is in the 
0.1% to 3.2% range, on average being equal to 1.3% for 
the (6-12) potential and 1.4% for the {6-exp) potential. 
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TABLE V: Comparison of experimental and calculated unit cell parameters for Fe(bpz)2(bipy) (at 1 atm). 





a, A 


6, A 


c, A 


A° 


V, 


angle," 


CM y/b 


H, kcal/mol 


The LS isomer 


calc. (6-12) 


16.319 


14.840 


10.685 


113.97 


2364 


92.55 


0.2699 


-49.78 


calc. (6-exp) 


16.136 


14.661 


10.697 


114.25 


2307 


91.93 


0.2724 


-40.22 


exp. (139 K) [42] 


16.086 


14.855 


10.812 


114.18 


2357 


90.91 


0.2754 




The HS isomer 


calc. (6-12) 


16.242 


15.178 


10.823 


113.60 


2445 


85.06 


0.2703 


-48.59 


calc. {6-exp) 


16.032 


14.995 


10.834 


113.92 


2381 


84.48 


0.2728 


-39.75 


exp. (293 K) [42] 


16.307 


15.075 


11.024 


114.95 


2457 


85.02 


0.2782 




exp. extrap. to 139 K 


16.16 


14.99 


11.04 


114.9 


2426-f-2429 








The difference between the HS and LS 


isomers 












calc. (6-12) 


-0.077 


0.338 


0.138 


-0.37 


81 


-7.49 


0.0004 


1.19 


calc. {6-exp) 


-0.104 


0.334 


0.137 


-0.33 


74 


-7.45 


0.0004 


0.47 


exp. extrap. to 139 K 


0.07 


0.14 


0.23 


0.7 


69-72 








exp. (30 K) [42 


-0.076 


0.347 


0.219 


1.09 


71.2 









As for the changes of these parameters in the course of the 
spin transition, in most cases the results predicted by the 
method of atom-atom potentials are in good agreement 
with the experimental data (the errors are ca. 0.05 — 0.1 
A). What is especially impressing is that the method is 
capable of reproducing decrease of some lattice periods 
in the course of the spin transition, which may happen in 
spite of the overall increase of the unit cell volume (the 
parameter b of the Fe(phen)2(NCS)2 crystal, the param- 
eter a of the Fe(bpz)2(bipy) crystal). However, we have 
three problematic cases: the variation of the parameter c 
of the Fe(phen)2(NCS)2 crystal (underestimated by the 
factor of 3-^5 times), and the variation of the parame- 
ters a and b of the Fe(btz)2(NCS)2 crystal (wrong sign 
of the result for a and underestimation by the factor of 
3-^4 times for 6). It is especially important that both 
Lennard-Jones and Buckingham potentials yield close re- 
sults. Trying to find an explanation for these errors, we 
noted that these three parameters are most sensitive to 
temperature changes. For example, the poorly predicted 
Ac of the Fe(phen)2(NCS)2 crystal (calculated from the c 
values extrapolated to the same temperature) changes in 
relative terms by 0.4% per 100 K, while both Aa and Ah 
- only by 0.2% per 100 K. Similarly, Aa, Ab and Ac of 
the Fe(btz)2(NCS)2 crystal decrease by 0.7%, 1.1% and 
0.6% per 100 K. This allows us to suggest that omission 
of the explicit treatment of the entropy contribution to 
the Gibbs energy, and thus the uncertainty in renormal- 
ization of the empirical parameters of the potentials, is 
one of the main sources of errors in the method in its 
current form, even if it is partially compensated by data 
correction for the thermal expansion. 

Another possible explanation (which docs not exclude 
the previous one) is that some specific interactions take 
place in these crystals, different from those in ordinary 
organic crystals used for fitting the presumably pure van 
der Waals interaction parameters. In this case, the per- 
formance of the method can be improved by correcting 
the parameters of atom-atom interactions. 




FIG. 8; Contacts in the crystal of Fe(plien)2(NCS)2. 



B. Contacts analysis and parameters adjustment 

To study this problem and yet further improve the 
performance of the method, we analyzed intermolecular 
contacts in the crystals, comparing atom-atom distances 
found in the experimental studies with those optimized 
with the parameters from [s^, [13] . The lists of the short- 
est atom-atom contacts (we selected those separated by 
less than the sum of the corresponding van der Waals 
radii) are given in Tables IVHIVIIII and they are also de- 
picted on Figs. [SlfTOl It is important to note that in all 
three materials the spin transition does not much affect 
the picture of intermolecular contacts. In other words, 
the shortest contacts in a LS crystal are also short (typi- 
cally, though not always, shorter than the sum of the van 
der Waals radii) contacts in its HS form, and vice versa. 

First of all, one can see that in most cases the shortest 
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TABLE VI: Shortest iritermolecular contacts in the LS and 
HS crystals of Fe(phen)2(NCS)2: interatomic distances (A) 
determined experimentally (exp.) or calculated with (6-12) 
and {6-exp) potentials and the sum of the van der Waals radii 
(vdW) of the atoms fs^. 



Pair 


R(exp.) 


R(vdW) 


R(6-12) 


R(6-ea;p) 


The LS isomer 


H-H 


2.093 


2.34 


2.374 


2.289 


S---C 


3.314 


3.55 


3.370 


3.341 


C---H 


2.589, 2.784 


2.92 


2.620, 2.801 


2.520, 2.866 


S--H 


2.832, 2.891, 


2.97 


3.162, 3.311, 


3.185, 3.419, 




2.911, 2.951 




3.052, 2.907 


2.983, 2.915 


The HS isomer 


H---H 


2.211 


2.34 


2.411 


2.307 


S---C 


3.357 


3.55 


3.345 


3.339 


C---H 


2.570, 2.750 


2.92 


2.635, 2.872 


2.509, 2.792 


S---H 


2.941 


2.97 


3.121 


3.156 



TABLE VH: Shortest intermolecular contacts in the LS and 
HS crystals of Fe(btz)2(NCS)2: interatomic distances (A) de- 
termined experimentally (exp.) or calculated with (6-12) and 
{6-exp) potentials and the sum of the van der Waals radii 
(vdW) of the atoms 



Pair 1 R(exp.) | R(vdW) | R(6-12) | R(6-ea;p) 


The LS isomer 


S---C 


3.275 


3.55 


3.352 


3.333 


S---H 


2.706, 2.743, 
2.942 


2.97 


2.886, 2.840, 
2.849 


2.755, 2.741, 
2.907 


C---H 


2.811 


2.92 


2.931 


2.805 


C---C 


3.453 


3.50 


3.668 


3.570 


The HS isomer 


S---C 


3.351 


3.55 


3.401 


3.384 


S---H 


2.893, 2.925 


2.97 


2.832, 2.904 


2.770, 2.805 


C---H 


2.848, 2.888 


2.92 


2.802, 2.851 


2.693, 2.755 



TABLE VIII: Shortest intermolecular contacts in the LS and 
HS crystals of Fe(bpz)2(bipy): interatomic distances (A) de- 
termined experimentally (exp.) or calculated with (6-12) and 
{6-exp) potentials and the sum of the van der Waals radii 
(vdW) of the atoms ;36]. 



Pair 


R(exp.) 


R(vdW) 


R(6-12) 


R(6-ea;p) 


The LS isomer 


C---H 


2.655, 2.658, 


2.92 


2.732, 2.716, 


2.651, 2.606, 




2.689, 2.817, 




2.815, 2.877, 


2.720, 2.763, 




2.879, 2.912 




2.895, 2.811 


2.835, 2.818 


H---H 


2.283, 2.332 


2.34 


2.387, 2.390 


2.301, 2.333 


C---C 


3.368, 3.374 


3.50 


3.327, 2.720 


3.387, 3.412 


The HS isomer 


C---H 


2.579, 2.719, 


2.92 


2.728, 2.739, 


2.622, 2.657, 




2.782, 2.813, 




3.002, 2.809, 


2.892, 2.707 




2.878, 2.900 




2.838, 3.189 


2.769, 3.040 


H-H 


2.318 


2.34 


2.477 


2.336 


C---C 


3.360, 3.420 


3.50 


3.250, 3.304 


3.194, 3.236 


N---H 


2.606 


2.67 


2.786 


2.642 




FIG. 9: Contacts in the crystal of Fe(btz)2(NCS)2. 




FIG. 10: Contacts in the crystal of Fe(bpz)2(bipy). 



contacts involve hydrogen atoms (S- - - H, C- - - H, N- - - H, 
or H- - - H). It is well known that coordinates of the hydro- 
gen atoms determined from X-ray diffraction may be sub- 
ject to significant errors unless tricks of crystallographic 
computing are used. While the X-H bond length is no- 
toriously underestimated due to the shift of the bonding 
electron pair towards the nonmetal X atom, the position 
of the X-H vector in three-dimensional space is correctly 
found. Thus, the H atom should "ride" on the nonmetal 
atom with a fixed bond length (e.g., C-H = 1.09 A, N-H 
= 1.01 A, 0-H = 0.96 A). Because the crystal structures 
under study seemingly did not profit from such "riding" 
H atoms approach, we may suggest that one of the main 
sources of mistakes in our results is the uncertainty in 
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the H positions. This also indicates that in the future 
research, when taking into consideration intramolecular 
degrees of freedom, one should take possible deformations 
of the C— H bonds into account. 

It is reasonable to suggest that the poorly described 
atom-atom interactions will be at the top of the list of the 
highest atom-atom repulsion energies. Indeed, if some in- 
teratomic distance increases when the system goes from 
the experimental configuration to the optimized one, 
the repulsion between the corresponding atoms weakens. 
Thus one can expect that the intensity of that interac- 
tion is overestimated, since such relaxation does not oc- 
cur in experiment. A similar reasoning applies to the 
strongest attractions as well. In practice, the picture is 
not so clear because molecules in organic crystals typi- 
cally have numerous contacts between various atoms. By 
analyzing the crystals formed by the Fe(phen)2(NCS)2 
or Fe(btz)2(NCS)2 molecules we found that the sulphur 
atoms play very important role in the intermolccular in- 
teractions (Fe(bpz)2(bipy) does not contain sulphur). As 
our calculations demonstrated, the S atoms participate in 
many close contacts with other atoms, thus providing a 
significant contribution to the repulsion within crystals; 
at the same time, their contributions to the attraction are 
also dominant (attraction energies of various S- ■ • S pairs 
are the largest by absolute value in these crystals; as for 
the S- ■ • C contacts, in some of them attraction is also 
very strong, while some other S- • ■ C contacts are among 
extreme cases of repulsion). 

The Fe(phen)2(NCS)2 molecule has S atoms only in 
the NCS groups while in the case of Fe(btz)2(NCS)2 the 
chelating ligand also contains the S atoms. We found that 
the S atoms of both types participate in the contacts with 
extremal values of the energy. Taking into consideration 
that the parameterization of the van der Waals energy of 
the S- • • X contacts (X — S, C, H) is not so well studied as 
compared to the C- ■ • C, C- • ■ H, and H- ■ • H interactions, 
and that some involvement of the lone pairs and vacant 
d-orbitals of the S atoms can complicate the approxima- 
tion of the S- • • X interactions by the center-symmetric 
atom-atom contributions, we suggest that improving the 
treatment of S- • • X (X = S, C, H) interaction energies 
may be another way of developing a better model of the 
atom-atom potentials for molecular crystals undergoing 
spin transitions. For example, the shortest C- ■ • S dis- 
tances are found to be ca. 0.2 A shorter than the sum of 
the van der Waals radii of the atoms. In the case of the 
S- • • H contacts, this contraction may reach even 0.27 A. 
Thus it is reasonable to suggest that due to some specific 
interactions, the optimal interatomic distances involving 
S atoms may be lower than determined by the standard 
parameterization. 

For that reason we adjusted the parameters for the 
S---C, S---S, S---H interactions in order to improve 
agreement between the experimental and modeled crys- 
tal lattice parameters. However, improvement of some 
of the calculated lattice parameters often increases the 
discrepancies for others. The situation is especially diffi- 



cult for the differences between the spin isomers Aa, Ab, 
Ac. Variation of the parameters for atom-atom contacts 
similarly affects the lattice parameters in the LS and HS 
crystals, hence the resulting change in those parameters 
is small and it can be only calculated rather than pre- 
dicted from any physical or geometrical reasoning. 

We performed a systematic quantitative study of the 
influence of the interaction parameters on the equilib- 
rium configurations of the crystals. To get the general 
understanding of this issue, we optimized the crystals of 
Fe(btz)2(NCS)2 with the interaction parameters slightly 
modified. We increased, one by one, parameters for each 
pair of atoms (the well depth and the equilibrium sep- 
aration) by 5% to estimate numerically the sensitivity 
of the energy contributions to the potential parameters. 
The choice of the Fe(btz)2(NCS)2 crystals was suggested 
by the fact that it is poorly described with the original 
parameterization: there are qualitative discrepancies for 
the changes in two out of three lattice parameters (Aa 
and Ab). Also we limited the consideration to the (6- 
exp) potential only, since the parameters of interaction 
between atoms of different elements were determined ex- 
plicitly [36j without any reference to the superposition 
approximation (except for the N- • ■ S pairs making little 
difference for the systems studied), and thus they can be 
varied separately. 

Changes in the optimal lattice parameters Sa, 6b, Sc, 
caused by 5% variations of each parameter of the atom- 
atom interaction energy, are given in Table HXl The table 
also specifies the differences between the experimental 
values of the lattice parameters and those calculated with 
the initial parameters of [3^ (exp.). As one can see from 
the numbers, most of the interaction parameters very 
slightly affect the optimal configuration of the crystals. 
Corrections caused by the well depth changes by 5% are 
a hundred times smaller than the difference between the 
experimental and calculated lattice parameters, leaving 
no hope to reduce the discrepancy by fitting the well 
depths within reasonable frames. The same applies to 
most of the atom-atom equilibrium separations on the 
corresponding interaction energy curves (C- • ■ C, H- ■ • H, 
N- • • N, N- • • H, etc.), though in this case the changes in 
the lattice parameters caused by a 5% increase of the 
distances are only tens times smaller than the required 
scale of correction. 

Only three parameters significantly affect the optimal 
structure of the crystal: the equilibrium separations for 
the S- • • C, S- • • H, and S- • • S pairs (hsted in the order of 
decreasing effect). This confirms our assumption made 
above on the basis of the interatomic contacts analysis 
that the contacts involving the S atoms need an improved 
treatment first of all. 

To do this, we performed a numerical minimization 
of the sum of squares of residuals / as a function of the 
equilibrium separations r for the S---C,S---H, and S- • ■ S 
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TABLE IX: Changes in the optimal lattice parameters a, b and c (10 '*A) of Fe(btz)2(NCS)2 caused by 5% increase of the 
{6-exp) potential parameters. 





The LS isomer 


The HS isomer 


The HS/LS difference 


Pair 


5a 1 Sb \ 5c 


5a \ 5b \ 5c 


5Aa\ 5Ab \ 5 Ac 


well depth 


H-H 


12 


14 


-3 


28 


-1 


-5 


16 


-15 


-2 


C-H 


-14 


-11 


1 


6 


-9 


11 


20 


2 


10 


N---H 


-13 


-5 


-13 


-15 


-6 


-13 


-2 


-1 





S---H 


126 


24 


6 


85 


41 


-10 


-41 


17 


-16 


C---C 


-28 


-20 


-19 


-27 


-20 


-14 


1 





5 


N---C 


-25 


-10 


-31 


-29 


-12 


-27 


-4 


-2 


4 


S---C 


-25 


27 


51 


-11 


18 


31 


14 


-9 


-20 


N---N 


-11 


-4 


-1 


-13 


-4 


3 


-2 





4 


S---S 


-23 


-13 


4 


-26 


-8 


19 


-3 


5 


15 


equilibrium distance 


H---H 


251 


216 


5 


483 


24 


-33 


232 


-192 


-38 


C---H 


396 


98 


365 


755 


74 


421 


359 


-24 


56 


N - -H 


-60 


-29 


-48 


-41 


-23 


-21 


19 


6 


27 


S---H 


2084 


828 


700 


1741 


991 


533 


-343 


163 


-167 


C---C 


168 


-55 


269 


277 


-61 


263 


109 


-6 


-6 


N---C 


-148 


-83 


-70 


-137 


-90 


-2 


11 


-7 


68 


S--C 


714 


1131 


2075 


894 


958 


1821 


180 


-173 


-254 


N---N 


-80 


-27 


271 


-109 


-18 


305 


-29 


9 


34 


S---S 


11 


-58 


1116 


-22 


28 


1226 


-33 


86 


110 


exp. 


705 


2930 


1354 


2098 


1843 


1131 


1393 


-1087 


-223 



atoms: 

/ = {a,LS,calc — aLS,expf + 

,calc ^LS.exp) 
iaHS,calc — aHS.exp) + {hnS^calc — bHS,exp) + (9) 
,calc CHS,exp ^calc ^^exp) H" 

{Abcalc - Abexp) +(A 

where Qcaic, bcaic, Ccaic stand for the optimal lattice pa- 
rameters calculated with the Buckingham potential pa- 
rameterization different from one in [3a | by rsci fSH and 
rss separations, and a^xp, bexp, Cexp are the experimen- 
tal lattice parameters (extrapolated to 293 K, if neces- 
sary). The result is that the equilibrium separation of 
the S- - ■ C contact should be increased by 6.6% (-1-0.26 
A), that of the S- - - H one - decreased by 2.3% (-0.08 
A), and that of the S- - - S one - decreased by 15% (—0.57 
A). Optimal values of the crystal lattice parameters of 
the Fe(btz)2(NCS)2 crystal, calculated with the adjusted 
parameters, are given in Table HVl fcalc. {6-exp) modif.). 

After fitting the S- - - C, S- - - H and S- - - S equilibrium 
distances all six unit cell dimensions became closer to the 
experimental values: the error in 6(LS) decreased from 
-0.32 A to -0.15 A, the error in a(HS) - from -0.21 
A to -0.10 A, in 6(HS) - from -0.18 A to -0.08 A, 
and so on. The values of the CM position and the ro- 
tation angle change, by contrast, insignificantly, in spite 
of the fact that they were not included in the treatment 
by the least squares method. The performance of the 
model in predicting the quantities Aa, Ab and Ac also 
significantly improved. The value of Aa shifted towards 



the experimental one and changed the sign to the correct 
one (positive instead of negative). A6 moved towards the 
experimental value, though this correction was only 1/3 
of the initial discrepancy. Finally, the Ac value shifted in 
the correct direction, but this time the change was even 
larger than the required one. An attempt to improve 
yet further the relation between the predicted and actual 
values of Aa, Ab, Ac and AV by another modification 
of the atom-atom interaction parameters (for example, 
by increasing the weights ascribed to the corresponding 
squares in the treatment by the least squares method) 
leads to catastrophic results for a, b and c of the pure 
LS and HS crystals: a tiny improvement by 0.01 A in 
Aa, Ab, Ac simultaneously leads to the growth in the 
discrepancies in a, b and c by ca. 0.1 A. 

We applied the same modified parameterization to the 
LS and HS crystals of Fe(phen)2(NCS)2. The results 
are given in Table iHll fcalc. {6-exp) modif.). In regard 
to the unit cell dimensions of the LS and HS crystals, 
the modification improved 4 out of 6 periods, especially 
those poorly described by the original parameterization: 
the error in 6(LS) decreased from -0.32 A to -0.19 A, 
in 6(HS) - from -0.29 A to -0.15 A. Error in unit cell 
volumes V decreased 2-3 times. At the same time, a sig- 
nificant error in the c value for the HS form appeared 
(0.33 A instead of 0.06 A). As for the changes in the lat- 
tice parameters, their values became more distant from 
the experimental values by 0.01 — 0.09 A. To sum up, the 
suggested modification generally improves the results of 
the model for both S-containing materials, though it fails 
to eliminate the errors completely. 
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TABLE X: Comparison of experimental and calculated unit 
cell parameters for Fe(plien)2(NCS)2 (at 1 GPa). 



System 




a, A 


b, A 


c, A 


V, P 


LS(1 GPa) 


calc. (6-12) 
calc. (^Q-exp) 
exp. (298 K) [48] 


13.060 
12.838 
12.656 


9.773 
9.700 
9.848 


17.183 
17.089 
16.597 


2193.2 
2128.0 
2068.6 


Difference, 
LS(1 GPa)/ 
HS(1 atm) 


calc. (6-12) 
calc. {6-exp) 
exp. 


-0.465 
-0.426 
-0.505 


-0.137 
-0.169 
-0.315 


-0.399 
-0.453 
-0.884 


-163.5 
-168.2 
-269.6 


Difference, 
LS(1 GPa)/ 
LS(1 atm) 


calc. (6-12) 
calc. {6-exp) 
exp. 


-0.125 
-0.155 
-0.114 


-0.149 
-0.161 
-0.242 


-0.163 
-0.193 
-0.625 


-76.0 
-86.0 
-150.5 



The adjustment of interaction parameters, described 
in this Subsection, does not claim to produce a new sys- 
tem of atom-atom parameters. We undertook it just to 
estimate how much improvement in the performance of 
the method at the expense of minor changes within the 
same theoretical paradigm may be done, and to illustrate 
that accurate treatment of the intermolecular contacts 
involving sulphur atoms are of primary importance for 
modeling the spin transition in S-containing materials. 



C. Contributions of intermolecular interactions to 
enthalpy 

The results described in the previous Sections demon- 
strate that the method of atom-atom potentials is capa- 
ble of modeling intermolecular interactions and reproduc- 
ing experimental data on the geometry of the unit cells. 
This allows us to go on to estimate the contributions of 
the van der Waals intermolecular forces to the energy 
(enthalpy) of the spin transitions, which cannot be ex- 
tracted from experimental data. The results are given in 
the last columns of Tables IIIHTVl First of all, one can see 
that this contribution may be either positive or neg ative, 
which corroborates the theoretical conclusion of dZ| . An- 
other important point is that the lattice contribution to 
the enthalpy of the spin transition is comparable with 
its total value. Though the estimates obtained with the 
Lennard-Jones and Buckingham potentials are somewhat 
different, the general picture is the same. For example, 
in the case of the Fe(phen)2(NCS)2 crystal we found this 
component to be equal to -1-1.81 kcal/mol (6-12) or -1-0.57 
kcal/mol {6-exp) or -fl.54 kcal/mol {6-exp modified), 
while the total experimental enthalpy (from the calori- 
metrical data) is -t-2.05 kcal/mol [9]. It means that one 
cannot neglect intermolecular interactions in calculating 
thermodynamical characteristics of the spin transitions 
in molecular crystals. (This conclusion was also made 
in [2^ on the basis of DFT calculations; however, the 
contribution of intermolecular interactions, which can be 
extracted from their results and ranging from 2 to 23 
kcal/mol, seems to be strongly overestimated). 



TABLE XL Comparison of experimental and calculated unit 
cell parameters for Fe(btz)2(NCS)2 (at 1 GPa). 



System 




a, A 


b, A 


c, A 


V, 


LS(1 GPa) 


calc. (6-12) 

C3,lc. (^Q-GXp^ 

exp. (298 K) [48] 


13.072 
12.877 
12.839 


10.410 
10.380 
10.454 


16.640 
16.502 
16.362 


2264.4 
2205.7 
2196.1 


Difference, 
LS(1 GPa)/ 
HS(1 atm) 


calc. (6-12) 
calc. {6-exp) 
exp. 


-0.171 
-0.366 
-0.449 


-0.313 
-0.343 
-0.407 


-0.307 
-0.445 
-0.558 


-142.2 
-200.9 
-245.8 


Difference, 
LS(1 GPa)/ 
LS(1 atm) 


calc. (6-12) 
calc. {6-exp) 
exp. 


-0.194 
-0.389 
-0.216 


-0.107 
-0.138 
-0.196 


-0.249 
-0.386 
-0.310 


-92.0 
-150.7 
-122.0 



D. Pressure effects 

Finally, we studied behavior of the crystal lattice 
parameters under the external hydrostatic pressure. 
Calculations were made for the Fe(phen)2(NCS)2 and 
Fe(btz)2(NCS)2 compounds, since the experimental data 
on the pressure effects on the spin transition are available 
only for these crystals [i^. We performed minimization 
of the enthalpies as a function of the lattice parameters, 
CM positions of the molecules, and their rotation an- 
gles, at two values of the pressure. The external pressure 
was accounted for by the PV term in the function to be 
minimized. The starting points of optimization were the 
experimental geometries. As previously, at the first step 
we minimized enthalpy as a function of five parameters, 
preserving the symmetry of the crystal, and after that 
checked that we get the global minima by allowing vari- 
ation of all 27 parameters mentioned above. The results 
for the lattice parameters of the LS forms at 1 GPa and 
298 K are given in Tables |Xl IXI) and the compressibil- 
ity coefficieirts at 1 atm and 1 GPa - in Tabic Kill As 
one can see from the tables, the high-pressure lattice pa- 
rameters are very well reproduced (errors are below 4%), 
though less accurately than those for the low pressure. 
As for the compressibility coefficients, in all the cases the 
correspondence between the calculated and experimental 
values is qualitative (the compressibility coefficients are 
underestimated by a factor of 1.5-^3 as compared to the 
experimental values). One can see that the Buckingham 
potential produces better values than the Lennard-Jones 
one. Taking into consideration that the {6-exp) parame- 
terization used in the present study is based only on the 
crystal structures measured at 1 atm, and very few con- 
tacts in those structures have distances corresponding to 
the repulsive branch of the potential (see Figs. 5-7 of Ref. 
^36^), we conclude that our results for the high-pressure 
structures are better than one could expect. 



V. CONCLUSION 

Numerical modeling of the spin transitions in molec- 
ular crystals is important from practical and theoretical 
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TABLE XIL Comparison of experimental and calculated com- 
pressibility coefficients (in 10"^ GPa"^) for Fe(phen)2(NCS)2 
and Fe(btz)2(NCS)2. 



System 






h 


kc 


kv 


Fe(phen)2(NCS)2 


calc. (6-12) 


0.14 


0.22 


0.14 


0.50 


HS, 1 atm 


calc. (6-exp) 


0.16 


0.22 


0.17 


0.56 




exp. (298 K) [48] 


0.21 


0.33 


0.53 


1.07 


Fefphen)ofNCS")-) 


calc. (6-12) 


0.07 


0.12 


0.07 


0.26 


LS, 1 GPa 


calc. i^-exp) 


0.09 


0.13 


0.09 


0.30 




exp. (298 K) [48] 


0.16 


0.28 


0.38 


0.82 


Fe(btz)2(NCS)2 


calc. (6-12) 


0.22 


0.14 


0.20 


0.56 


HS, 1 atm 


calc. (6-exp) 


0.25 


0.14 


0.18 


0.57 




exp. (298 K) [48] 


0.41 


0.43 


0.37 


1.21 


Fe(btz)2(NCS)2 


calc. (6-12) 


0.11 


0.08 


0.11 


0.30 


LS, 1 GPa 


calc. i^-exp) 


0.13 


0.09 


0.11 


0.34 




exp. (298 K) [48] 


0.28 


0.33 


0.28 


0.89 



viewpoints. There is no alternative to calculations explic- 
itly taking into account the composition and structure 
of interacting molecules (instead of representing them 
by spheres, or ellipsoids, or octahedra etc., immersed in 
an elastic media), both for the purposes of theoretical 
study of the transition mechanisms and for prediction 
of phenomenological parameters for macroscopic models. 
Meanwhile, the modern quantum chemical methods are 
hardly applicable to such objects, because their accuracy 
level is not sufficient to calculate the required values (for 
example, enthalpies of the spin transitions). 

We demonstrate that the atom-atom potentials can be 
used for analysis of intermolecular contributions to the 
structure and energy of spin-active crystals. Indeed, in- 
termolecular contacts in these crystals are those between 
the C, H, N, S, etc. organogenic atoms, while the metal 
atom and its bonds with the donor atoms of the ligands 
are hidden inside the complex. As a consequence of that, 
the van der Waals interactions in the spin-active crystals 
can be approximated similarly to those in ordinary or- 
ganic molecular crystals. In the present paper we checked 
for the first time the possibility to use the atom-atom po- 
tentials method for this class of objects. 

In all the cases the symmetry groups of optimized 
crystals coincided with those found in experiment; the 
unit cell volumes were calculated with the precision of 
0.5 — 4%. Errors in the predicted lattice parameters did 
not exceed 3% at the ambient pressure and 4% at 1 GPa. 
Direction (sign) and magnitude of the changes of the lat- 
tice parameters and molecules positions in the unit cell in 



the course of the temperature- and pressure-driven spin 
transitions were reproduced correctly. The compress- 
ibility coefficients are in a qualitative agreement with 
their experimental values, although 1.5-^3 times underes- 
timated. Thus the accuracy of the method of atom-atom 
potentials is quite sufficient at the present level of the 
theory. We attempted to improve the parameterization, 
which is based on the {6-exp) parameters from [36] and 
differs from it in the S- ■ • H, S- • • C and S- • ■ S equilibrium 
separations. The results of this fitting of the parameters 
demonstrate that the performance of the method can be 
significantly improved by adjustment to the specific cases 
under study, and that the energy of interactions involving 
sulphur atoms is the crucial term for adequate treatment 
of spin transitions in the crystals studied. 

Our study shows that any reliable calculation of spin 
transition parameters (such as transition enthalpy) must 
take into account intermolecular interactions. According 
to our estimates for the Fe(phen)2(NCS)2 crystal, the 
van der Waals contribution to the transition enthalpy is 
about -1-0.6-^1.8 kcal/mol (as compared with the total 
transition enthalpy of -1-2.05 kcal/mol). 

We believe that the accuracy of the method used in this 
paper is limited by (i) implicit treatment of the entropy 
effects (through fitting the interaction parameters, rather 
than explicit calculation of frequencies of intermolecular 
oscillations); (ii) uncertainty of the H atoms positions 
in the experimental X-ray structures; (iii) description of 
the energy of interactions involving the S atoms (due to 
possible involvement of lone pairs and vacant d-orbitals 
of the sulphur atoms). Nevertheless, even the current 
level of precision is enough for using the method of atom- 
atom potentials to study the spin transitions in molecular 
crystals. 
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